###################
# CREATE BASELINE #
###################

# Running Time: ~ 30 minutes

rm(list=ls())
starts <- 5
KK = 4221


######################################################
# Import Demand Estimation Results and Export to csv #
######################################################

for (k in 1:KK) { 
objective_MPEC_starts = matrix(0,1,starts)
status_MPEC_starts = matrix(0,1,starts)  
filename <- sprintf("2_results_full/resultsr_%s_%s.RData",k ,1) 
if (file.exists(filename)) {
load(filename)
X_MPEC_starts = X_MPEC 
objective_MPEC_starts[1,1] = objective_MPEC
status_MPEC_starts[1,1] = status_MPEC 
for (reps in 2:starts) {  
  filename <- sprintf("2_results_full/resultsr_%s_%s.RData",k ,reps )
  load(filename)
  X_MPEC_starts <- cbind(X_MPEC_starts, X_MPEC)
  objective_MPEC_starts[1,reps] = objective_MPEC
  status_MPEC_starts[1,reps] = status_MPEC 
}
X_MPEC = X_MPEC_starts
objective_MPEC = objective_MPEC_starts
status_MPEC = status_MPEC_starts
filename_X_MPEC = sprintf("3_baseline/X_MPEC_%s.csv", k)
write.csv(X_MPEC, file = filename_X_MPEC)
filename_objective_MPEC = sprintf("3_baseline/objective_MPEC_%s.csv", k)
write.csv(objective_MPEC, file = filename_objective_MPEC)
filename_status_MPEC = sprintf("3_baseline/status_MPEC_%s.csv", k)
write.csv(status_MPEC, file = filename_status_MPEC)
filename_data = sprintf("3_baseline/data_%s.csv", k)
write.csv(data, file = filename_data)
filename_FE = sprintf("3_baseline/FE_%s.csv", k)
write.csv(FE, file = filename_FE)
filename_v_save = sprintf("3_baseline/v_save_%s.csv", k)
write.csv(v_save, file = filename_v_save)
filename_Y_save = sprintf("3_baseline/Y_save_%s.csv", k)
write.csv(Y_save, file = filename_Y_save)
}
}

write.csv(K, file = "3_baseline/K.csv")                                      # Save K to File
write.csv(K_A, file = "3_baseline/K_A.csv")                                  # Save K_A to File


load("2_results_full/income_mu_sigma.RData")                 # Import Demographics
names(income_mu_sigma) <- c("declarant", "year", "mu", "sigma")
write.csv(income_mu_sigma, file = "3_baseline/income_mu_sigma.csv", row.names=F)        # Export Demographics to csv

######################################################################################################################


#####################################
# Determine Products that converged #
#####################################

converged <- matrix(0,KK,1)  
for (k in 1:KK) { 
  status_MPEC_starts = matrix(0,1,starts)  
  filename <- sprintf("2_results_full/resultsr_%s_%s.RData",k ,1)  
  if (file.exists(filename)) {
  load(filename)
  status_MPEC_starts[1,1] = status_MPEC 
  for (reps in 2:starts) {  
    filename <- sprintf("2_results_full/resultsr_%s_%s.RData",k ,reps )
    load(filename)
    status_MPEC_starts[1,reps] = status_MPEC 
  }
  status_MPEC = status_MPEC_starts
  converged[k,1] = any(status_MPEC == 0)                                                                  # Test if any iteration converged (status == 0)
  }
}

#########################################################################################################################################################


#####################################
# Replace results if converged == 0 #
#####################################

converged_final <- converged

for (k in 1:KK) { 
  if (converged[k,1] == 0) {
  objective_MPEC_starts = matrix(0,1,starts)
  status_MPEC_starts = matrix(0,1,starts)  
  filename <- sprintf("1_results_subsample/resultsr_%s_%s.RData",k ,1)  
  if (file.exists(filename)) {
  load(filename)
  X_MPEC_starts = X_MPEC 
  objective_MPEC_starts[1,1] = objective_MPEC
  status_MPEC_starts[1,1] = status_MPEC 
  for (reps in 2:starts) {  
    filename <- sprintf("1_results_subsample/resultsr_%s_%s.RData",k ,reps)  
    load(filename)
    X_MPEC_starts <- cbind(X_MPEC_starts, X_MPEC)
    objective_MPEC_starts[1,reps] = objective_MPEC
    status_MPEC_starts[1,reps] = status_MPEC 
  }
  X_MPEC = X_MPEC_starts
  objective_MPEC = objective_MPEC_starts
  status_MPEC = status_MPEC_starts
  filename_X_MPEC = sprintf("3_baseline/X_MPEC_%s.csv", k)
  write.csv(X_MPEC, file = filename_X_MPEC)
  filename_objective_MPEC = sprintf("3_baseline/objective_MPEC_%s.csv", k)
  write.csv(objective_MPEC, file = filename_objective_MPEC)
  filename_status_MPEC = sprintf("3_baseline/status_MPEC_%s.csv", k)
  write.csv(status_MPEC, file = filename_status_MPEC)
  filename_data = sprintf("3_baseline/data_%s.csv", k)
  write.csv(data, file = filename_data)
  filename_FE = sprintf("3_baseline/FE_%s.csv", k)
  write.csv(FE, file = filename_FE)
  filename_v_save = sprintf("3_baseline/v_save_%s.csv", k)
  write.csv(v_save, file = filename_v_save)
  filename_Y_save = sprintf("3_baseline/Y_save_%s.csv", k)
  write.csv(Y_save, file = filename_Y_save)
  converged_final[k,1] <- any(status_MPEC == 0) 
  }
  }
}



pc8plus_i_k <- matrix(0, KK, 2)

i_k <- 0
for (k in 1:KK) { 
  if (converged_final[k,1] == 1) {
    i_k = i_k + 1
    pc8plus_i_k[k,1] = k
    pc8plus_i_k[k,2] = i_k
    filename_X_MPEC = sprintf("3_baseline/X_MPEC_%s.csv", k)
    filename_X_MPEC2 = sprintf("3_baseline/X_MPEC_%s.csv", i_k)
    file.rename(filename_X_MPEC, filename_X_MPEC2)
    filename_objective_MPEC = sprintf("3_baseline/objective_MPEC_%s.csv", k)
    filename_objective_MPEC2 = sprintf("3_baseline/objective_MPEC_%s.csv", i_k)
    file.rename(filename_objective_MPEC, filename_objective_MPEC2)
    filename_status_MPEC = sprintf("3_baseline/status_MPEC_%s.csv", k)
    filename_status_MPEC2 = sprintf("3_baseline/status_MPEC_%s.csv", i_k)
    file.rename(filename_status_MPEC, filename_status_MPEC2)
    filename_data = sprintf("3_baseline/data_%s.csv", k)
    filename_data2 = sprintf("3_baseline/data_%s.csv", i_k)
    file.rename(filename_data, filename_data2)
    filename_FE = sprintf("3_baseline/FE_%s.csv", k)
    filename_FE2 = sprintf("3_baseline/FE_%s.csv", i_k)
    file.rename(filename_FE, filename_FE2)
    filename_v_save = sprintf("3_baseline/v_save_%s.csv", k)
    filename_v_save2 = sprintf("3_baseline/v_save_%s.csv", i_k)
    file.rename(filename_v_save, filename_v_save2)
    filename_Y_save = sprintf("3_baseline/Y_save_%s.csv", k)
    filename_Y_save2 = sprintf("3_baseline/Y_save_%s.csv", i_k)
    file.rename(filename_Y_save, filename_Y_save2)
  }
}
  

write.csv(file = 'pc8plus_i_k.csv', pc8plus_i_k)


non_converged <- c(3473, 3488, 3504, 3506, 3509, 3512, 3646, 3766, 3768, 3811, 3832, 3833, 3844, 3870:3876, 3881:3882, 3885, 3888, 4021, 4023, 4025, 4027:4030, 4032, 4132:4147, 4166:4170, 4190:4199)

for (i in non_converged) {
filename_X_MPEC = sprintf("3_baseline/X_MPEC_%s.csv", i)
file.remove(filename_X_MPEC) 
filename_objective_MPEC = sprintf("3_baseline/objective_MPEC_%s.csv", i)
file.remove(filename_objective_MPEC)
filename_status_MPEC = sprintf("3_baseline/status_MPEC_%s.csv", i)
file.remove(filename_status_MPEC)
filename_data = sprintf("3_baseline/data_%s.csv", i)
file.remove(filename_data)
filename_FE = sprintf("3_baseline/FE_%s.csv", i)
file.remove(filename_FE)
filename_v_save = sprintf("3_baseline/v_save_%s.csv", i)
file.remove(filename_v_save)
filename_Y_save = sprintf("3_baseline/Y_save_%s.csv", i)
file.remove(filename_Y_save)
}

